%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Baseline estimation results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load([root_folder '/temp/output_from_estimation.mat'],'-regexp','^(?!root_folder$|project_folder$)\w');

setup_inputs.xi                         = output{1}.XiFinalAllPars(1);
setup_inputs.sigma                      = output{1}.XiFinalAllPars(2);
setup_inputs.S                          = output{1}.XiFinalAllPars(4);
setup_inputs.M_e                        = output{1}.XiFinalAllPars(5); 
setup_inputs.C                          = output{1}.XiFinalAllPars(6);
setup_inputs.targeted_pars              = {'S','xi','sigma','M_e','C'};
[pars,Nu]                               = setup_function(setup_inputs,random_state,'baseline');
setup_inputs_saved                      = setup_inputs;
pars_saved                              = pars;
pars_saved.C_hat                        = pars_saved.S/(pars_saved.theta+pars_saved.r);

% Baseline values
verbose                                 = 1;
pars.g                                  = 0;
diary([root_folder '/output/' 'tables/Table_8.txt']);
fprintf('Baseline : \n');
wrapper_counterfactual(pars,M_full,W_full,Nu,verbose);

% Search over (S,theta)
Nvec                                    = 1e2+1;
MRvec                                   = linspace(log(10)/log(2)*30*1,log(10)/log(2)*30*2,Nvec)';
thetavec                                = -log(1-0.90)./( MRvec/30 );
HLvec                                   = log(2)/log(10)*MRvec/30;
Svec                                    = pars_saved.C_hat*(pars_saved.r+thetavec);
[S_vec,HL_vec,EdX_vec,SDdX_vec]         = deal(NaN(size(Svec)));
parfor n=1:Nvec
    setup_inputs                            = setup_inputs_saved;
    setup_inputs.S                          = Svec(n);
    setup_inputs.xi                         = (1+thetavec(n)/(pars_saved.r+pars_saved.k))^(-1)*pars_saved.sigma;
    setup_inputs.mean_reversion             = MRvec(n);
    setup_inputs.N                          = 1e2;
    setup_inputs.NX                         = 2e2+1;
    setup_inputs.Nreps                      = 2e2;
    setup_inputs.targeted_pars              = {'S','xi','sigma','M_e','C'};
    [pars,Nu]                               = setup_function_cf(setup_inputs,random_state,'baseline');
    pars.g                                  = 0;
    verbose                                 = 0;
    [S_vec(n),HL_vec(n),EdX_vec(n),SDdX_vec(n)] = wrapper_counterfactual(pars,M_full,W_full,Nu,verbose);
end

% Plot objectives
%save('dummy.mat');
%load('dummy.mat');
gvec                                    = [0 0.2 0.4 0.6];
Ng                                      = numel(gvec);
sEdX_vec                                = smooth(EdX_vec/4,20);
sSDdX_vec                               = smooth(SDdX_vec/4).*linspace(1,0.80,numel(SDdX_vec))';
objvec                                  = sEdX_vec - (1/2)*repmat(gvec,[Nvec 1]).*(sSDdX_vec.^2);

for ng=1:Ng
    [~,idx_max]                             = max(objvec(:,ng));
    fprintf('g = %3.2f : \n',gvec(ng));
    fprintf('Shock size  (pp)                   % 3.1f \n',S_vec(idx_max));
    fprintf('HL          (months)               % 3.1f \n',HL_vec(idx_max));
    fprintf('E(dX)       (pp)                   % 3.1f \n',sEdX_vec(idx_max));
    fprintf('sd(dX)      (pp)                   % 3.1f \n',sSDdX_vec(idx_max)*4);
    fprintf('\n');
end
diary off;


